library(rdd)
library(rdpower)
library(dummies)
library(Matching)
library(Rmisc)
library(extrafont)
library(magrittr)
library(dplyr)
#install.packages('extrafontdb')
#system.file("fontmap", "fonttable.csv", package="extrafontdb")
#font_import(paths = NULL, recursive = TRUE, prompt = TRUE,pattern = NULL)
#loadfonts()
library(ggplot2)
library(stargazer)
###Replicating stata robust cluster standard errors (Jens Hainmueller code)
require(sandwich)
require(lmtest)
library(plm)
library(MASS)
library(quantmod)
options(error=NULL)
library(dummies)
options(warn=2, error=recover)
#options(warn=0)
library(rdrobust)
set.seed(1234)

#Purpose of file: loading necessary data, packages, and defining functions for replication for "Anticipating Dissent: The Repression of Politicians in Pinochet's Chile"

#R version 3.5.3 (2019-03-11) macOS Mojave 10.14.6

##Note: bootstrapping takes several hours to run


##robust clustered standard errors (replicates stata, code by Jens Hainmueller)
vcovCluster <- function(model, cluster){
  require(sandwich)
  require(lmtest)
  if(nrow(model.matrix(model))!=length(cluster)){
    stop("check your data: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)          
  K <- model$rank  
  if(M<50){
    warning("Fewer than 50 clusters, variances may be unreliable (could try block bootstrap instead).")
  }
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}


####multiple plots in the same panel
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  require(grid)
  plots <- c(list(...), plotlist)
  numPlots = length(plots)
  if (is.null(layout)) {
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    for (i in 1:numPlots) {
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

setwd("/Users/janeesberg/Dropbox/PhD/RepressedPols/RepressedPols_May2019")

###1973 elections
dips73<-read.csv('dips73_cleaned_JOP.csv')

##target list
contreras<-read.csv('contreras_list_coded.csv')

##1969 elections
dips69<-read.csv('dips69_cleaned_JOP.csv')

##regidores
regidores_full<-read.csv('regidores_cleaned_JOP.csv')

#################################################################
##calculating quotient_diff for d'Hondt proportional system######
##################################################################

#rank ordering politicians by party and district
rankorder<-function(df){
  df<-df[order(-df$Votos),]
  df$rank<-1:length(df$list)
  return(df)
}

dips73<-ddply(dips73, .(list, Distrito), rankorder)

##check it correctly ordered
cbind(dips73$rank[dips73$Distrito=="Rancagua, Caupolican, San Vicente"], 
      dips73$list[dips73$Distrito=="Rancagua, Caupolican, San Vicente"],
      dips73$Votos[dips73$Distrito=="Rancagua, Caupolican, San Vicente"])

##get quotients
dips73$quotient<-(dips73$vote_total/dips73$rank)

dips73<-ddply(dips73, .(Distrito), transform, quotient_perc=quotient/sum(quotient)*100)
dips73<-ddply(dips73, .(Distrito), transform, quotient_diff_winners=quotient_perc-max(quotient_perc[elected==0]))
dips73<-ddply(dips73, .(Distrito), transform, quotient_diff_losers=quotient_perc-min(quotient_perc[elected==1]))

dips73$quotient_diff<-ifelse(dips73$quotient_diff_winners>0, dips73$quotient_diff_winners, dips73$quotient_diff_losers)
dips73$quotient_diff<-dips73$quotient_diff/2

table(I(dips73$quotient_diff>0)[dips73$elected==1])
table(I(dips73$quotient_diff<0)[dips73$elected==0])



##Defining different forms of dataset
dips73$UP<-NA
dips73$UP[dips73$list=="UP"]<-1
dips73$UP[is.na(dips73$UP)]<-0

dips73_nousopo<-subset(dips73, list!="USOPO")
dips73_left<-subset(dips73, list!="CODE")
dips73$leftDC<-ifelse(dips73$list=="UP", 1, 0)
dips73$leftDC<-ifelse(dips73$Partido=="PDC", 1, dips73$leftDC)
dips73$leftDC<-ifelse(dips73$Partido=="USOPO", 1, dips73$leftDC)
dips73_leftDC<-subset(dips73, dips73$leftDC==1)
dips73_leftDC_nousopo<-subset(dips73_leftDC, dips73_leftDC$list!="USOPO")
dips73_noasylum<-subset(dips73, CD!="ASYLUM_EXILE")
dips73_nousopo_noasylum<-subset(dips73_nousopo, CD!="ASYLUM_EXILE")
dips73_leftDC_nousopo_noasylum<-subset(dips73_leftDC_nousopo, CD!="ASYLUM_EXILE")



###produces estimates for elected and unelected candidates and difference
ttest_electedvunelected<-function(type, dataset){
  ttest_90<-t.test(type[dataset$elected==1], type[dataset$elected==0], conf.level=.9)
  elected<-ttest_90$estimate[1]*100
  unelected<-ttest_90$estimate[2]*100
  difference<-(ttest_90$estimate[1]-ttest_90$estimate[2])*100
  lower_90<-ttest_90$conf.int[1]*100
  upper_90<-ttest_90$conf.int[2]*100
  
  ttest_95<-t.test(type[dataset$elected==1], type[dataset$elected==0], conf.level=.95)
  lower_95<-ttest_95$conf.int[1]*100
  upper_95<-ttest_95$conf.int[2]*100
  
  data.frame(rbind(cbind(elected, unelected, difference, lower_90,
                         upper_90, '90%'), cbind(elected, unelected, difference, lower_95, upper_95, '95%')))
}


####function to bootstrap logit coefficients
dips73$agrupacion_factor<-as.factor(dips73$agrupacion)
dips73$list_factor<-as.factor(dips73$list)

dips73_nousopo$agrupacion_factor<-as.factor(dips73_nousopo$agrupacion)
dips73_nousopo$list_factor<-as.factor(dips73_nousopo$list)

dips73_leftDC_nousopo$agrupacion_factor<-as.factor(dips73_leftDC_nousopo$agrupacion)
dips73_leftDC_nousopo$list_factor<-as.factor(dips73_leftDC_nousopo$list)

coef=list()
logitCoefficients<-function(df, variables, iterations){
  district<-unique(df$district)
  n=length(district)
  
  dataset<-cbind(df[variables], df['district'])
  
  for (j in 1:iterations){
    v=sample(length(district), n, replace=TRUE)
    newdata=NULL
    
    for (i in 1:n){
      newdata=rbind(newdata, subset(dataset, district==v[i]))
    }
    
    if ('anyrepression' %in% colnames(newdata)){
      model<-tryCatch(glm(anyrepression~., data=subset(newdata, select=-c(district)), family='binomial'),error=function(e) NULL)
      if(is.null(model)) {next}
      else{coef[[j]]=coef(model)[1:2]}
    }
    else if ('victim' %in% colnames(newdata)){
      model<-tryCatch(glm(victim~., data=subset(newdata, select=-c(district)), family='binomial'),error=function(e) NULL)
      if(is.null(model)) {next}
      else{coef[[j]]=coef(model)[1:2]}
    }
    else if ('search' %in% colnames(newdata)){
      model<-tryCatch(glm(search~., data=subset(newdata, select=-c(district)), family='binomial'),error=function(e) NULL)
      if(is.null(model)) {next}
      else{coef[[j]]=coef(model)[1:2]}
    }
  }
  coef
}


####Getting marginal effects using observed values approach for logistic regression
marginalEffectsCalc<-function(model, data, coefficient, modelparts){
  data$constant<-model$coefficients[1]
  data$iv<-coefficient
  effect0<-exp(data$constant+data$modelparts)/(1+exp(data$constant+data$modelparts))
  effect1<-exp(data$constant+data$iv+data$modelparts)/(1+exp(data$constant+data$iv+data$modelparts))
  effect<-effect1-effect0
  mean(effect)
}


###OLS coefficients

OLS_Coefficients<-function(df, variables, iterations){
  district<-unique(df$district)
  n=length(district)
  
  dataset<-cbind(df[variables], df['district'])
  
  for (j in 1:iterations){
    v=sample(length(district), n, replace=TRUE)
    newdata=NULL
    
    for (i in 1:n){
      newdata=rbind(newdata, subset(dataset, district==v[i]))
    }
    
    if ('anyrepression' %in% colnames(newdata)){
      model<-tryCatch(lm(anyrepression~., data=subset(newdata, select=-c(district))),error=function(e) NULL)
      if(is.null(model)) {next}
      else{coef[[j]]=coef(model)[1:2]}
    }
    else if ('victim' %in% colnames(newdata)){
      model<-tryCatch(lm(victim~., data=subset(newdata, select=-c(district))),error=function(e) NULL)
      if(is.null(model)) {next}
      else{coef[[j]]=coef(model)[1:2]}
    }
    else if ('search' %in% colnames(newdata)){
      model<-tryCatch(lm(search~., data=subset(newdata, select=-c(district))),error=function(e) NULL)
      if(is.null(model)) {next}
      else{coef[[j]]=coef(model)[1:2]}
    }
  }
  coef
}


###For other DVs
logitCoefficientsAlternative<-function(dataset, iterations, string, var){
  district<-unique(dataset$district)
  dataset$iv<-unlist(dataset[var])
  districtdummies<-dummy(dataset$agrupacion)
  colnames(districtdummies)<-gsub(string, 'agrupacion', colnames(districtdummies))
  dataset<-cbind(dataset, districtdummies)
  
  n=length(district)
  iterations=iterations
  
  for (j in 1:iterations){
    v=sample(length(district), n, replace=TRUE)
    newdata=NULL
    
    for (i in 1:n){
      rownames(newdata) <- c()
      rownames(dataset) <- c()
      newdata=rbind(newdata, subset(dataset, dataset$district==v[i]))
    }
    
    error1<-tryCatch(glm(victim~iv+factor(UP)+agrupacion2+agrupacion3+
                           agrupacion4+agrupacion5+agrupacion6+agrupacion7+
                           agrupacion8+agrupacion9+agrupacion10, data=newdata, family="binomial"), error=function(e) e)
    if (class(error1)[1]=='simpleError'){
      next
    }
    else{
      rownames(newdata) <- c()
      reg1=glm(victim~iv, data=newdata, family='binomial')
      reg1p=glm(victim~iv+factor(UP), data=newdata, family="binomial")
      reg1pd=glm(victim~iv+factor(UP)+agrupacion2+agrupacion3+
                   agrupacion4+agrupacion5+agrupacion6+agrupacion7+
                   agrupacion8+agrupacion9+agrupacion10, data=newdata, family="binomial")
      alt1_coef[[j]]=coef(reg1)
      alt1p_coef[[j]]=coef(reg1p)
      alt1pd_coef[[j]]=coef(reg1pd)
      
    }
  }
  c(alt1_coef, alt1p_coef, alt1pd_coef)
}


####forest plots
dips73$levels<-ifelse(dips73$quotient_diff<quantile(dips73$quotient_diff[dips73$elected==0], .5), 0, 1)
dips73$levels<-ifelse(dips73$quotient_diff>quantile(dips73$quotient_diff[dips73$elected==1], .5), 3, dips73$levels)
dips73$levels<-ifelse(dips73$elected==1&dips73$levels==1, 2, dips73$levels)

dips73_nousopo$levels<-ifelse(dips73_nousopo$quotient_diff<quantile(dips73_nousopo$quotient_diff[dips73_nousopo$elected==0], .5), 0, 1)
dips73_nousopo$levels<-ifelse(dips73_nousopo$quotient_diff>quantile(dips73_nousopo$quotient_diff[dips73_nousopo$elected==1], .5), 3, dips73_nousopo$levels)
dips73_nousopo$levels<-ifelse(dips73_nousopo$elected==1&dips73_nousopo$levels==1, 2, dips73_nousopo$levels)

dips73_leftDC_nousopo$levels<-ifelse(dips73_leftDC_nousopo$quotient_diff<quantile(dips73_leftDC_nousopo$quotient_diff[dips73_leftDC_nousopo$elected==0], .5), 0, 1)
dips73_leftDC_nousopo$levels<-ifelse(dips73_leftDC_nousopo$quotient_diff>quantile(dips73_leftDC_nousopo$quotient_diff[dips73_leftDC_nousopo$elected==1], .5), 3, dips73_leftDC_nousopo$levels)
dips73_leftDC_nousopo$levels<-ifelse(dips73_leftDC_nousopo$elected==1&dips73_leftDC_nousopo$levels==1, 2, dips73_leftDC_nousopo$levels)

whiskerPlotLevels<-function(dataframe, dv, xlab, ylab, ggtitle){
  ttest_level0<-t.test(dv[dataframe$levels==0], conf.level=.9)
  mean_level0<-ttest_level0$est
  lower_level0<-ttest_level0$conf.int[1]
  upper_level0<-ttest_level0$conf.int[2]
  
  ttest_level0_95<-t.test(dv[dataframe$levels==0], conf.level=.95)
  mean_level0_95<-ttest_level0_95$est
  lower_level0_95<-ttest_level0_95$conf.int[1]
  upper_level0_95<-ttest_level0_95$conf.int[2]
  
  ttest_level1<-t.test(dv[dataframe$levels==1], conf.level=.9)
  mean_level1<-ttest_level1$est
  lower_level1<-ttest_level1$conf.int[1]
  upper_level1<-ttest_level1$conf.int[2]
  
  ttest_level1_95<-t.test(dv[dataframe$levels==1], conf.level=.95)
  mean_level1_95<-ttest_level1_95$est
  lower_level1_95<-ttest_level1_95$conf.int[1]
  upper_level1_95<-ttest_level1_95$conf.int[2]
  
  ttest_level2<-t.test(dv[dataframe$levels==2], conf.level=.9)
  mean_level2<-ttest_level2$est
  lower_level2<-ttest_level2$conf.int[1]
  upper_level2<-ttest_level2$conf.int[2]
  
  ttest_level2_95<-t.test(dv[dataframe$levels==2], conf.level=.95)
  mean_level2_95<-ttest_level2_95$est
  lower_level2_95<-ttest_level2_95$conf.int[1]
  upper_level2_95<-ttest_level2_95$conf.int[2]
  
  ttest_level3<-t.test(dv[dataframe$levels==3], conf.level=.9)
  mean_level3<-ttest_level3$est
  lower_level3<-ttest_level3$conf.int[1]
  upper_level3<-ttest_level3$conf.int[2]
  
  ttest_level3_95<-t.test(dv[dataframe$levels==3], conf.level=.95)
  mean_level3_95<-ttest_level3_95$est
  lower_level3_95<-ttest_level3_95$conf.int[1]
  upper_level3_95<-ttest_level3_95$conf.int[2]
  
  mean<-rbind(mean_level0, mean_level1, mean_level2, mean_level3)
  lower<-rbind(lower_level0, lower_level1, lower_level2, lower_level3)
  upper<-rbind(upper_level0, upper_level1, upper_level2, upper_level3)
  lower95<-rbind(lower_level0_95, lower_level1_95, lower_level2_95, lower_level3_95)
  upper95<-rbind(upper_level0_95, upper_level1_95, upper_level2_95, upper_level3_95)
  
  level=rbind(1, 2, 3, 4)
  df<-data.frame(cbind(level, mean, lower, upper))
  df$CI<-'90%'
  df95<-data.frame(cbind(level, mean, lower95, upper95))
  df95$CI<-'95%'
  df<-rbind(df, df95)
  colnames(df)<-c("level", "mean", "lower", "upper", 'CI')
  
  ggplot(data=df, aes(x=factor(level), y=mean, ymin=lower, ymax=upper, col=factor(level), linetype=CI)) +
    geom_pointrange(size=.75) + 
    xlab(xlab) + ylab(ylab) +
    theme_bw() +  theme(legend.position="none", axis.title=element_text(size=20),
                        text=element_text(family="Times New Roman", size=20),
                        plot.title = element_text(size=20),
                        panel.border=element_rect(colour='black', fill=NA, size=.75))+
    scale_color_manual(values=c("dark grey", "black", "black", "dark grey"))+
    scale_linetype_manual(values=c("solid", "dotted"))+
    ggtitle(ggtitle)+
    ylim(0, .7)+
    theme(plot.title = element_text(hjust = 0.5))
}


###bootstrap local linear regression
LLR=NA
LLRhalf=NA
LLR2=NA
LLRparty=NA
LLRpartyhalf=NA
LLRparty2=NA

bootstrapLLR<-function(data, iterations){
  set.seed(822)
  distritos_num<-unique(data$cluster)
  n=length(distritos_num)
  for (j in 1:iterations){
    v=sample(distritos_num, n, replace=TRUE)
    newdata=NULL
    for (i in 1:n){
      newdata=rbind(newdata, subset(data, data$cluster==v[i]))
    }
    mod1<-lm(Y ~ Tr + Xl + Xr, weights = w, data = subset(newdata, w > 0))
    mod2<-lm(Y ~ Tr + Xl + Xr, weights=whalf, data=subset(newdata, whalf>0))
    mod3<-lm(Y ~ Tr + Xl + Xr, weights=w2, data=subset(newdata, w2>0))
    
    mod4<-lm(Y ~ Tr + Xl + Xr+PartidoDR+PartidoIC+
               PartidoMAPU+PartidoPADENA+PartidoPCCh+PartidoPDC+PartidoPIR+PartidoPN+PartidoPR+
               PartidoPS, weights = w, data = subset(newdata, w > 0))
    mod5<-lm(Y ~ Tr + Xl + Xr+PartidoDR+PartidoIC+
               PartidoMAPU+PartidoPADENA+PartidoPCCh+PartidoPDC+PartidoPIR+PartidoPN+PartidoPR+
               PartidoPS, weights = whalf, data = subset(newdata, whalf > 0))
    mod6<-lm(Y ~ Tr + Xl + Xr+PartidoDR+PartidoIC+
               PartidoMAPU+PartidoPADENA+PartidoPCCh+PartidoPDC+PartidoPIR+PartidoPN+PartidoPR+
               PartidoPS, weights = w2, data = subset(newdata, w2 > 0))
    LLR[j]=mod1$coefficients[2]
    LLRhalf[j]=mod2$coefficients[2]
    LLR2[j]=mod3$coefficients[2]
    LLRparty[j]=mod4$coefficients[2]
    LLRpartyhalf[j]=mod5$coefficients[2]
    LLRparty2[j]=mod6$coefficients[2]
    
  }
  data.frame(LLR, LLRhalf, LLR2, LLRparty, LLRpartyhalf, LLRparty2)
}

###prepare data for LLR
makeLLR<-function(originalData, dependentvariable, independentvariable, cutpoint, string){
  originalData$distritos_num<-as.numeric(originalData$Distrito)
  X<-independentvariable-cutpoint
  Xl <- (X < 0) * X
  Xr <- (X >= 0) * X
  Tr <- as.integer(X >= 0)
  Y<-dependentvariable
  cluster<-originalData$distritos_num
  bw=IKbandwidth(X = X, Y = Y, cutpoint = 0, kernel = "triangular")
  w<-kernelwts(X, 0, bw=bw, kernel="triangular")
  whalf<-kernelwts(X, 0, bw=bw/2, kernel="triangular")
  w2<-kernelwts(X, 0, bw=2*bw, kernel="triangular")
  dummies<-data.frame(dummy(originalData$Partido))
  colnames(dummies)<-gsub(string, 'Partido', colnames(dummies))
  PartidoDR<-dummies$PartidoDR
  PartidoIC<-dummies$PartidoIC
  PartidoMAPU<-dummies$PartidoMAPU
  PartidoPADENA<-dummies$PartidoPADENA
  PartidoPCCh<-dummies$PartidoPCCh
  PartidoPDC<-dummies$PartidoPDC
  PartidoPIR<-dummies$PartidoPIR
  PartidoPN<-dummies$PartidoPN
  PartidoPR<-dummies$PartidoPR
  PartidoPS<-dummies$PartidoPS
  cluster <- originalData$distritos_num
  data.frame(Y, Tr, X, Xl, Xr, w, whalf, w2, cluster,PartidoDR, PartidoIC,PartidoMAPU, PartidoPADENA, PartidoPCCh, PartidoPDC, PartidoPIR, PartidoPN,PartidoPR, PartidoPS)
}
makeAlt<-function(originalData, dependentvariable, independentvariable, cutpoint, string, bw){
  originalData$distritos_num<-as.numeric(originalData$Distrito)
  X<-independentvariable-cutpoint
  Xl <- (X < 0) * X
  Xr <- (X >= 0) * X
  Tr <- as.integer(X >= 0)
  Y<-dependentvariable
  cluster<-originalData$distritos_num
  w<-kernelwts(X, 0, bw=bw, kernel="triangular")
  whalf<-kernelwts(X, 0, bw=bw/2, kernel="triangular")
  w2<-kernelwts(X, 0, bw=2*bw, kernel="triangular")
  dummies<-data.frame(dummy(originalData$Partido))
  colnames(dummies)<-gsub(string, 'Partido', colnames(dummies))
  PartidoDR<-dummies$PartidoDR
  PartidoIC<-dummies$PartidoIC
  PartidoMAPU<-dummies$PartidoMAPU
  PartidoPADENA<-dummies$PartidoPADENA
  PartidoPCCh<-dummies$PartidoPCCh
  PartidoPDC<-dummies$PartidoPDC
  PartidoPIR<-dummies$PartidoPIR
  PartidoPN<-dummies$PartidoPN
  PartidoPR<-dummies$PartidoPR
  PartidoPS<-dummies$PartidoPS
  cluster <- originalData$distritos_num
  data.frame(Y, Tr, X, Xl, Xr, w, whalf, w2, cluster,PartidoDR, PartidoIC,PartidoMAPU, PartidoPADENA, PartidoPCCh, PartidoPDC, PartidoPIR, PartidoPN,PartidoPR, PartidoPS)
}

####Testing for balance
covariateBalanceTestData<-function(originalData, dependentvariable){
  originalData$distritos_num<-as.numeric(originalData$Distrito)
  X <-originalData$quotient_diff
  Xl <- (X < 0) * X
  Xr <- (X >= 0) * X
  Tr <- as.integer(X >= 0)
  Y <- dependentvariable
  cluster<-originalData$distritos_num
  bw=IKbandwidth(X = X, Y = Y, cutpoint = 0, kernel = "triangular")
  w<-kernelwts(originalData$quotient_diff, 0, bw=bw, kernel="triangular")
  data.frame(Y, Tr, X, Xl, Xr, w, cluster)
}
covariate_late=NA
covariateBalance<-function(data){
  distritos_num<-unique(data$cluster)
  n=length(distritos_num)
  
  for (j in 1:iterations){
    v=sample(distritos_num, n, replace=TRUE)
    newdata=NULL
    
    for (i in 1:n){
      newdata=rbind(newdata, subset(data, data$cluster==v[i]))
    }
    mod<-lm(Y ~ Tr + Xl + Xr, weights = w, data = subset(newdata, w > 0))
    
    covariate_late[j]=mod$coefficients[2]
  }
  covariate_late
}

##linear regression RDD
ols_anyrepression=NA
ols_anyrepression_party=NA
ols_physicalcoercion=NA
ols_physicalcoercion_party=NA

linearRegressionFun<-function(data, trim, string){
  data$distritos_num<-as.numeric(data$Distrito)
  distritos_num<-unique(data$distritos_num)
  partidodummies<-dummy(data$Partido)
  colnames(partidodummies)<-gsub(string, 'Partido', colnames(partidodummies))
  data<-cbind(data, partidodummies)
  n<-length(distritos_num)
  
  trim<-subset(data, abs(quotient_diff)<trim)
  for (j in 1:iterations){
    v=sample(length(distritos_num), n, replace=TRUE)
    newdata=NULL
    
    for (i in 1:n){
      newdata=rbind(newdata, subset(trim, data$distritos_num==v[i]))
    }
    reg_anyrepression=lm(anyrepression~factor(elected)*quotient_diff, data=newdata)
    reg_anyrepression_party=lm(anyrepression~factor(elected)*quotient_diff+
                                 PartidoDR+PartidoIC+PartidoMAPU+PartidoPADENA+PartidoPCCh+
                                 PartidoPDC+PartidoPIR+PartidoPN+PartidoPR+
                                 PartidoPS+PartidoUSOPO, data=newdata)
    reg_physicalcoercion=lm(victim~factor(elected)*quotient_diff, data=newdata)
    reg_physicalcoercion_party=lm(victim~factor(elected)*quotient_diff+
                                    PartidoDR+PartidoIC+PartidoMAPU+PartidoPADENA+PartidoPCCh+
                                    PartidoPDC+PartidoPIR+PartidoPN+PartidoPR+
                                    PartidoPS+PartidoUSOPO, data=newdata)
    ols_anyrepression[j]=reg_anyrepression$coefficients[2]
    ols_anyrepression_party[j]=reg_anyrepression_party$coefficients[2]
    ols_physicalcoercion[j]=reg_physicalcoercion$coefficients[2]
    ols_physicalcoercion_party[j]=reg_physicalcoercion_party$coefficients[2]
  }
  data.frame(ols_anyrepression, ols_anyrepression_party, ols_physicalcoercion, ols_physicalcoercion_party)
}


####1969 functions

length(dips69$agrupacion)==length(dips69$Distrito)

rankfunction2<-function(df){
  df<-df[order(-df$Votos),]
  df$rank<-1:length(df$Partido)
  return(df)
}

dips69<-ddply(dips69, .(Partido, Distrito), rankfunction2)
dips69$quotient<-dips69$vote_total/dips69$rank
dips69<-ddply(dips69, .(Distrito), transform, quotient_perc=quotient/sum(quotient)*100)
dips69<-ddply(dips69, .(Distrito), transform, quotient_diff_winners=quotient_perc-max(quotient_perc[elected==0]))
dips69<-ddply(dips69, .(Distrito), transform, quotient_diff_losers=quotient_perc-min(quotient_perc[elected==1]))

dips69$quotient_diff<-ifelse(dips69$quotient_diff_winners>0, dips69$quotient_diff_winners, dips69$quotient_diff_losers)
dips69$quotient_diff<-dips69$quotient_diff/2

dips69_leftDC<-subset(dips69, Partido!="DN")
dips69_leftDC<-subset(dips69_leftDC, Partido!="PN")
dips69_leftDC<-subset(dips69_leftDC, Partido!="IN")

dips69_nounelected<-subset(dips69, Partido!="IN")
dips69_nounelected<-subset(dips69_nounelected, Partido!="DN")
dips69_nounelected<-subset(dips69_nounelected, Partido!="USP")
dips69_nounelected<-subset(dips69_nounelected, Partido!="PSD")

dips69_leftDC_nounelected<-subset(dips69_nounelected, Partido!="PN")



